/*
 * EquilibriumSolver.cpp
 *
 *  Created on: 26 Jul 2011
 *      Author: Allan
 */

#include "EquilibriumSolver.h"

// C++ includes
#include <cmath>
#include <functional>
#include <iomanip>
#include <memory>
using namespace std;
using namespace std::placeholders;

// Boost includes
#include <boost/foreach.hpp>
using namespace boost;

// GeoReact includes
#include "Utils/Algorithms.h"
#include "Utils/SpeciesMap.h"
#include "Utils/SpeciesUtils.h"

EquilibriumSolver::EquilibriumSolver() :
tolerance(1.0E-6), maxiterations(100)
{}

EquilibriumSolver::EquilibriumSolver(const Multiphase& multiphase, const ReactionSystem& reactions) :
tolerance(1.0E-6), maxiterations(100)
{
	Initialiaze(multiphase, reactions);
}

void EquilibriumSolver::Initialiaze(const Multiphase& multiphase, const ReactionSystem& eReactions)
{
	// Set data members
	this->multiphase = multiphase;
	this->eReactions = eReactions;
	
	// Initialiaze the equilibrium problem instance
	problem.Initialiaze(multiphase, eReactions);
	
	// Initialiaze the multiroot solver (note that the modified Newton-Raphson algorithm works too, but the algorithms Hybridj and Hybridsj does not)
	solver.Initialiaze(multiphase.GetNumSpecies(), StandardNewton);
	solver.SetSolverOptions(tolerance, maxiterations);
}

void EquilibriumSolver::SetEquilibriumConditions(const vector<EquilibriumCondition>& eConditions)
{
	problem.SetEquilibriumConditions(eConditions);
}

void EquilibriumSolver::SetSolverOptions(double tolerance, unsigned maxiterations)
{
	// Set data members
	this->tolerance = tolerance;
	this->maxiterations = maxiterations;
	
	solver.SetSolverOptions(tolerance, maxiterations);
}

const EquilibriumSolver::Diagnostic EquilibriumSolver::Solve(double T, double P, VectorXd& n)
{
	// Set the paramenters of the equilibrium problem
	problem.SetParameters(T, P);
	
	// Set the diagnostic instance
	Diagnostic diagnostic = {1, false};
	
	// The solution vector (n) at the previous iteration (used for application of under relaxation to the solution)
	VectorXd nOld = n;
	
	// Solve the equilibrium problem
	while(solver.Iterate(problem, n) && diagnostic.iterations < maxiterations)
	{
		// Apply the under relaxation to the vector (n) to avoid negative values of compositions
		UnderRelax(n, nOld);
		
		// Update the vector (nOld)
		nOld = n;
		
		// Update (iteration)
		++diagnostic.iterations;
	}
	
	// Check convergence of the previous iterations
	diagnostic.converged = diagnostic.iterations < maxiterations;
	
	return diagnostic;
}

void EquilibriumSolver::UnderRelax(VectorXd& n, const VectorXd& nOld) const
{
	for(unsigned i = 0; i < unsigned(n.rows()); ++i)
	{
		const double delta  = n[i] - nOld[i];
		const double factor = max(1.0, -2.0 * delta/nOld[i]);
		
		n[i] = nOld[i] + delta/factor;
	}
}

